Granular clustering in a hydrodynamic simulation 
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We present a numerical simulation of a granular material using hydrodynamic equations. We 
show that, in the absence of external forces, such a system phase-separates into high density and 
low density regions. We show that this separation is dependent on the inelasticity of collisions, 
and comment on the mechanism for this clustering behavior. Our results are compatible with the 
granular clustering seen in experiments and molecular dynamic simulations of inelastic hard disks. 



One of the key differences between a granular material 
and a regular fluid is that the grains of the former lose 
energy with each collision, while the molecules of the lat- 
ter do not. Even when the inelasticity of the collisions is 
small, it can give rise to dramatic effects, including the 
Maxwell Demon effect|Q and, the topic of this paper, the 
phenomenon of granular clustering. Experiments^, |j| 
and molecular dynamics simulations Q] alike show that 
granular gases in the absence of gravity do not become 
homogeneous with time, but instead form dense clusters 
of stationary particles surrounded by a lower density re- 
gion of more energetic particles. From a particulate point 
of view, one can explain these clusters by noting that 
when a particle enters a region of slightly higher den- 
sity, it has more collisions, loses more energy, and so is 
less able to leave that region, thus increasing the local 
density and making it more likely for the next particle 
to be captured. We, however, are interested in describ- 
ing this behavior using hydrodynamics. There is con- 
siderable workjs) deriving granular hydrodynamics from 
kinetic theory, focusing on analytical treatments of the 
long-wavelength behavior of the system. Goldhirsch and 
ZannettiQ, for instance, describe clustering as the result 
of a hydrodynamic instability: a region of slightly higher 
density has more collisions, and thus has a lower tem- 
perature. A lower temperature means a lower pressure, 
which attracts more mass from the surrounding higher- 
pressure region. Their long-wavelength stability analysis 
shows that in a system of hydrodynamic equations simi- 
lar to Eq. [j], higher-density regions do indeed have lower 
pressure, fueling the instability. Our approach here is 
different. In this paper, we discuss whether a coarse- 
grained description, in terms of local particle, momen- 
tum, and energy densities, can be used to treat charac- 
teristic behaviors of granular materials as a self-contained 
dynamical system. The spirit is similar to the use of den- 
sity functional theory in the treatment of freezing j6|. We 
present the results of numerical simulations for such a 
description in the case of clustering, and show that when 
inelasticity is present, one indeed sees the system phase- 
separate into regions of high and low density. 

We begin with a number density field p, a flow velocity 
field u, and a "temperature" field T. These are related 
by a standard set of hydrodynamic equations for granular 



materials [|8|: 

dp 

dt 
d(puj) 

dt 
dT 

~di~ 



-Vi(pui) 

-VjP - Vj(puiUj) + V 'j(r)ijkiV 'kUi) 

-Vi(u,-T) + -Vj(KVjT) 
P 

+~f?yfcl(V 'iUj)(V \u{) - jT. (1) 



where repeated indices are summed over, and where P 
is the pressure, k is the bare thermal conductivity, and 
rji-jki = rj{8ik5 3 i +S t i5kj +Sij6ki) is the isotropic bare vis- 
cosity tensor. These equations bear much in common 
with those for normal fluids^]. The most important ad- 
dition is that of the inelastic term, — 7T. The parameter 
7 is a measure of the inelasticity of the collisions, and 
is proportional to (1 — e 2 ), where e is the coefficient of 
restitution. Using kinetic theory results the trans- 
port coefficients are chosen to depend on temperature 
and density: 



K = k T 1/2 
r, = ^pT 1 ' 2 
7 = 7oT 1/2 . 



(2) 



Typically, work in granular hydrodynamics is done 
in low-density regimes, where grains may be treated as 
point particles interacting via collisions. When simulat- 
ing aggregation, however, one must take excluded volume 
into account. We do this by introducing a barrier in the 
pressure P(p) at some maximum (close-packed) density 
Pq. This is in addition to the usual hydrodynamic pres- 
sure pT. We choose in particular the simple quadratic 
form 



P = P T+U(p 2 -pl)6(p-p ), 



(3) 



where U is a positive parameter, 9{x) is the unit step 
function, and po is the close-packed density. We intro- 
duced this method in an earlier paper and it is a 
simple way to model the incompressibility of the system 
at high densities Q. 

We evaluate our equations in two dimensions using a 
finite-difference Runge-Kutta method, on a 100 x 100 
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square lattice with periodic boundary conditions. (See 
footnote |Q for more details.) The lattice spacing is cho- 
sen to be large enough so that each site contains a number 
of grains, and we can consider the density to be a contin- 
uous variable. We start with random initial conditions 
p = 0.1 + 0.001ri(z, x), u z — r%{z, x), u x — r%(z, x), and 
T = 1 + 0.1ri(z,x), where Ti{z,x) are random numbers 
chosen between —1 and 1. The model's other parame- 
ters for the data presented here are 770 = 50, Ko = 1, 
U = 4 x 10 4 , po = 0.2, and 70 = 50 or 0. All values are 
in arbitrary units, with a time step of At = 10 -3 . 
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FIG. 2: The width of the density distributions, simply de- 
fined as the difference between the maximum density and the 
minimum density, for 70 = and 70 = 50. 
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FIG. 1: The density distribution at t = 800 when 70 = 50, 
clearly showing phase separation. The top picture is a spatial 
distribution of the density; darker shades correspond to higher 
densities. The lower graph is a histogram of the density, with 
bin size A(p/po) = 0.005. The initial distribution is shown 
by the light gray bar in the center. 

Figure |l| shows the density distribution for the system 
after it has evolved to a time t = 800. One can clearly 
see that, from a narrow initial distribution centered at 
p/ Po = 0.5, the system has separated into regions of high 



and low density, with relatively sharp borders. Calculat- 
ing the structure factor S(r) = (Sp(r')Sp(r' + r)) , shows 
that the density is correlated out to 8 or 9 lattice spac- 
ings. To show that this clustering behavior is due to 
the inelastic term in Eq. ^, we duplicate the run with 
70 = 0. Figure || shows the width of the density distri- 
bution as a function of time for the inelastic and elastic 
cases. Without inelasticity, the system quickly becomes 
homogeneous. 

This clustering behavior of the system appears to be 
robust. The parameter values chosen for this run were 
not optimized, and a variety of alternate choices result in 
the same qualitative effect. The time over which cluster- 
ing occurs does depend on different parameters, and in 
cases where the viscosity is too large, the density distribu- 
tion does not form the peaks at its extremes, but instead 
remains broad and relatively flat. This also occurs when 
one removes the pressure barrier at pQ. In both cases, 
however, there is distinct and permanent phase separa- 
tion. 

It is clear that our simulations reproduce clustering. 
As mentioned in the introduction, it is believed 0] that 
clustering occurs because of an anticorrelation between 
density and temperature: regions of high density have 
more collisions, and thus more dissipation, and thus less 
total energy. Fig. || shows the density distributions cor- 
responding to the highest and the lowest temperatures. 
Notice that the high-temperature distribution is missing 
the high-density peak present in the other two graphs. 
On the other hand, in the low-temperature distribution, 
the high-density peak is relatively larger than in the gen- 
eral case. Both points suggest that high-density regions 
tend to be at low temperatures, which agrees with the 
common wisdom. However, the anticorrelation is not as 
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FIG. 4: The spatial distribution of the log of the temperature 
at t — 800. Just as before, darker shades correspond to higher 
temperatures. Compare this figure with Fig. [j] to look for 
correlations between density and temperature. 



FIG. 3: Density distributions for the 70 = 50 case at t = 800. 
The middle graph is the distribution for the entire system. 
The upper graph is the distribution only for those sites on 
the lattice where the temperature is less than 5 x 10~ 7 . (The 
average temperature at this time is 1.3 x 10~ 6 .) Similarly, the 
lower graph is the distribution for those sites with tempera- 
ture greater than 5 x 10 -6 . 

clear-cut as one might expect: there is no correlation be- 
tween low densities and high temperatures, for instance, 
and even the high density sites have a broad range of tem- 
peratures which only tend toward the low side. Fig. || is 
a density plot of the temperature. When one compares 
this with Fig. [|, there appears to be more structure to 
the system than the simple explanation would suggest. 

One interesting point about the system is that it takes 
a long time to come to a stop. Fig. || shows the distri- 
bution of density at four different times. Also shown is 
the average kinetic energy of the system over time; it is 
decaying at a slower-than-exponential rate. This occurs 
because the viscosity, being proportional to the tempera- 
ture, is decaying to zero, and is thus very small by the end 
of this run. It is difficult to say whether this feature cor- 
responds to the results of experiments and other simula- 
tions, because the nature of freezing in gravity- free granu- 
lar systems is generally subtle. In event-driven molecular 
dynamic simulations, one finds the phenomenon of inelas- 
tic collapse pl| , where particles on the verge of freezing in 



a cluster suffer an infinite number of collisions in a finite 
amount of time. In experiment^], there is the opposite 
effect: one must vibrate real two-dimensional systems of 
particles to create clustering, as the surface friction will 
otherwise freeze the system before any clustering can oc- 
cur. If desired, one could bring our model to a halt with 
the introduction of a surface friction term, or by replac- 
ing the current temperature dependence of the viscosity 
with 

7?c<77o ( o(T + T min ) 1 / 2 . (4) 

In summary, we have recreated in hydrodynamic simu- 
lation the clustering behavior in granular materials which 
was predicted by hydrodynamic theory and observed in 
kinetic simulations. We show that this behavior is di- 
rectly dependent on the inelasticity parameter of our 
equations. This supports the notion that granular mate- 
rials can be, to some extent at least, described by means 
of coarse-grained variables in a nonlinear hydrodynamical 
setting. However, we also suggest that long-wavelength 
hydrodynamics alone may not be sufficient to fully de- 
scribe granular behavior, such as the actual mechanism 
for granular clustering. 

We thank Professor Todd DuPont for his assistance 
in constructing our simulations, and Professors Hcinrich 
Jaeger, Sidney Nagel, and Thomas Witten for helpful 
conversations. This work was supported by the Materials 
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FIG. 5: At the top, images of the spatial density distribution 
at four times, using the same key as in Fig. hi Note that 
the clusters continue to move around for long periods of time. 
Below these, a plot showing that the average kinetic energy 
(Ekin = (p( u z+ w x))) is decaying at a slower-than-exponential 
rate. 
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